Partial suppression of the radial orbit instability in stellar systems 
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ABSTRACT 

It is well known that the simple criterion proposed originally by Polyachenko 
and Shukhman (1981) for the onset of the radial orbit instability, although being 
generally a useful tool, faces significant exceptions both on the side of mildly 
anisotropic systems (with some that can be proved to be unstable) and on the 
side of strongly anisotropic models (with some that can be shown to be stable). 
In this paper we address two issues: Are there processes of collisionless collapse 
that can lead to equilibria of the exceptional type? What is the intrinsic struc- 
tural property that is responsible for the sometimes noted exceptional stability 
behavior? To clarify these issues, we have performed a series of simulations of 
collisionless collapse that start from homogeneous, highly symmetrized, cold ini- 
tial conditions and, because of such special conditions, are characterized by very 
little mixing. For these runs, the end-states can be associated with large values 
of the global pressure anisotropy parameter up to 2Kr/KT ~ 2.75. The highly 
anisotropic equilibrium states thus constructed show no significant traces of ra- 
dial anisotropy in their central region, with a very sharp transition to a radially 
anisotropic envelope occurring well inside the half- mass radius (around 0.2 tm)- 
To check whether the existence of such almost perfectly isotropic "nucleus" might 
be responsible for the apparent suppression of the radial orbit instability, we 
could not resort to equilibrium models with the above characteristics and with 



^present address: Space Telescope Science Institute, 3700 San Martin Drive Baltimore MD 21218 USA 



- 2 - 



analytically available distribution function; instead, we studied and confirmed 
the stability of configurations with those characteristics by initializing N-body 
approximate equilibria (with given density and pressure anisotropy profiles) with 
the help of the Jeans equations. 

Subject headings: stellar dynamics — galaxies: evolution — galaxies: kinematics 
and dynamics — galaxies: structure — methods: n-body simulations 

1. Introduction 

The study of the stability of a system of which the distribution function is available 
in analytical form can be carried out along three different approaches: N-body simulations 
(e.g., see Henon 1973; Merritt & Aguilar 1985; Barnes et al. 1986; Aguilar & Merritt 1990; 
Allen et al. 1990; Stiavelli & Sparke 1991; and following papers), hnear modal analysis 
(Polyachenko & Shukhman 1981; Palmer & Papaloizou 1987, 1988; Weinberg 1989, 1991, 
1994; Saha 1991, 1992; Bertin et al. 1994; see also Pridman & Polyachenko 1984; Palmer 
1993, and references therein), and energy principles (e.g., Sygnet et al. 1984; Kandrup & 
Sygnet 1985; Goodman 1988, and references therein). Spherical stellar systems with too 
many radial orbits have thus been found to be unstable to perturbations that break the 
spherical symmetry and remove the excess of kinetic energy in the radial degree of freedom 
(see Fridman & Polyachenko 1984 and Palmer 1993). The existence of such radial orbit 
instability was first investigated by Polyachenko & Shukhman (1981), who proposed an 
"empirical criterion" for the onset of instability {2Kr/KT > 1.7 ± 0.25) based on the global 
content of kinetic energy in the radial with respect to that in the tangential degrees of 
freedom. Following investigations (Barnes 1985; Merritt & Aguilar 1985; Aguilar & Merritt 
1990; Palmer 1993) noted that different families of models may exhibit different anisotropy 
thresholds for the instability, thus widening the uncertainty interval around the value of 1.7 
suggested earlier. For example, while for the so-called /oo and f^'^^ models a threshold value 
similar to that suggested by Polyachenko and Shukhman may be applicable (see Bertin & 
StiaveUi 1989 and Trenti & Bertin 2005, hereafter TB05), for the family of Dehnen (1993) 
density profiles with an Osipkov- Merritt type (Osipkov 1979; Merritt 1985) of anisotropy 
profile it has been argued that such value is as high as 2.5 (Meza & Zamorano 1997). 
On the other hand, systems with an arbitrarily small content of radial anisotropy can be 
unstable, although with very small growth rates (as shown by Palmer & Papaloizou (1987) 
by means of a linear modal analysis). 

The radial orbit instability is thought to play an important role during the formation 
of self-gravitating structures from coUisionless collapse via incomplete violent relaxation, a 



-3- 



process that has been argued to be the main driver for the formation of eUiptical galaxies (van 
Albada 1982). In fact, if a system starts from sufficiently cold initial conditions (i.e. with 
a low initial virial ratio u = {2K/\W\)t=o ^ 0.15), it collapses with stars falling in almost 
radially toward the center, often ending up in triaxial configurations, the origin of which is 
attributed to the radial orbit instability (Palmer et al. 1990; Udry 1993; Hjorth & Madsen 
1995). In general, when such coUisionless collapse leads to "realistic" final configurations, 
the values of pressure anisotropy achieved at the end of the collapse turn out to be consistent 
with the threshold value associated with the instability criterion proposed by Polyachenko 
and Shukhman (1981). Currently, galaxy formation is approached in the generally accepted 
cosmological context of hierarchical clustering (see e.g. Meza et al. 2003, and references 
therein), but the mechanism of incomplete violent relaxation remains an important ingredient 
and thus quantifying the effects of the radial orbit instability is relevant to the explanation of 
the observed fiattening of gravitational structures in cosmological simulations (see Hopkins 
et al. 2005). 

In this paper we present the results of some numerical experiments that are aimed 
at clarifying the following two issues: Are there processes of coUisionless collapse able to 
lead to equilibria of the exceptional type, that is equilibria violating the criterion proposed 
by Polyachenko and Shukhman (1981)? What is the structural property that makes some 
models with relatively low levels of global radial anisotropy unstable and others with high 
levels of global radial anisotropy stable? 

As to the first issue, in a set of simulations of coUisionless collapse (Trenti et al. 2005, 
hereafter TBvA05) we have investigated the role of mixing in phase space on the relaxation 
of the end products of the simulations. We have thus run, for comparison, a number of ex- 
periments in which phase mixing is inefficient. These are simulations with highly symmetric 
initial conditions. In spite of the cold initial conditions, with u < 0.1, and of the high level 
of radial anisotropy achieved in the final configurations, with 2Kj./Kt sometimes above 2.5 
(during the process of collapse values of 2Kj./Kt above 10 are reached), the radial orbit 
instability did not develop. From inspection of the structure of the final products of these 
simulations, it appears that most likely the physical factor at the basis of the observed par- 
tial suppression of the radial orbit instability is the existence of an almost perfectly isotropic 
central region that is realized in these systems, which appears to be more efficient for more 
concentrated models. This is the clue that we consider in order to address the second issue 
raised above. 

In fact, we note that at fixed content of global anisotropy 2Kt./Kt-, the f'^"^ and f^o 
models (for which the threshold value of 2Kr/KT is close to 1.7) are less isotropic locally, in 
their core, than similar systems with anisotropy profiles of the type introduced by Osipkov 
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(1979) and Merritt (1985), as studied by Meza & Zamorano (1997), which in turn are less 
isotropic in their central regions than the equilibrium states obtained at the end of our set 
of simulations. In the opposite direction, we recall that the generalized polytropic spheres, 
with small content of global anisotropy, found to be unstable by Palmer & Papaloizou (1987) 
are characterized by a constant and finite local anisotropy level down to r = 0. To confirm 
this picture and to study the effects on the radial orbit instability of central density and 
anisotropy profiles decoupled from each other, we cannot resort to equilibrium models with 
the above characteristics and with analytically available distribution function. We have 
then constructed with the help of the Jeans equations N-body equilibria with density and 
anisotropy profiles qualitatively similar to those found as end states for the simulations of 
collisionless collapse violating the criterion proposed by Polyachenko & Shukhman (1981) 
and thus demonstrated that some configurations with 2Kj./Kt up to ~ 2.9 do not show 
evidence of rapid evolution. 

The paper is organized as follows. In the next Section wc describe some special processes 
of collisionless collapse leading to highly anisotropic quasi-equilibrium states. In Sect. 3 we 
investigate the evolution of candidate collisionless equilibrium configurations constructed 
numerically starting from the Jeans equations and identify a quasi-equilibrium state with 
2Kr/ Kt ~ 2.9. Conclusions are presented in Sect. 4. In the Appendix we outline the method 
used to generate approximate equilibrium configurations with given density and anisotropy 
profiles starting from the Jeans equations. 

2. Exceptionally stable equilibria resulting from special processes of 

collisionless collapse 

Most of the numerical simulations discussed in this paper have been carried out with the 
recently developed particle-mesh code described in TBvAOS (for further tests and details see 
also Trenti 2005). The code solves the Poisson equation by expanding the density and the 
potential in spherical harmonics and is well suited for the study of the stability of collisionless 
systems (TB05). We have also run a few simulations with the fast tree code GyrFalcON 
(Dehnen 2000, 2002) to check, with success, that our results are not biased by the specific 
choice of the numerical code. 

As briefiy noted in the Introduction, it is well known (see van Albada 1982 and many 
following papers) that the collapse of a cold cloud of stars can lead, on a dynamical time- 
scale, to the formation of an equilibrium structure characterized by a quasi-isotropic core 
and a radially anisotropic halo with projected density profiles similar to the i?^/^ law (de 
Vaucouleurs 1948). The mechanism of incomplete violent relaxation (Lynden-Bell 1967) has 
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thus been advocated as a possible scenario for the formation of eUiptical galaxies. Several 
investigations (van Albada 1982; McGlynn 1984; Londrillo et al. 1991) have shown that the 
outcome of the simulations generally resembles observed systems provided the initial condi- 
tions are clumpy. This condition finds a natural interpretation in the current cosmological 
framework for structure formation, as the use of clumpy initial conditions in collisionlcss col- 
lapse simulations may effectively mimic some aspects of hierarchical merging. Simulations 
of mergers of stellar systems, i.e. of incomplete violent relaxation, are generally considered 
as an important tool to study the properties of elliptical galaxies, not only in view of an 
explanation of the origin of the universality of the observed surface brightness profiles, but 
also in relation to the origin and evolution of the Fundamental Plane (Nipoti et al. 2003; 
Gonzalez-Garci'a & van Albada 2003). Observational evidence is now accumulating in sup- 
port of the importance of mergers as a formation mechanism for elliptical galaxies (e.g. see 
van Dokkum & Ellis 2003; Bell et al. 2004), although many arguments and new findings 
tend to point to a merger process in which dissipation plays a significant role. In addition, 
Khochfar & Burkert (2003) (see also Jesseit et al. 2005) argue from semi-analytic modeling 
that 50 % of present days ellipticals are relics of mergers between spheroidal systems. 

For the present study, we focus on a particular set of simulations that start from homoge- 
neous spherically symmetric configurations obtained by symmetrizing a clumpy configuration 
with ten cold clumps, the kinetic energy of which is in the collective motion of their center 
of mass. In most of the runs u = {2K/\W\)t=o < 0.1; the clump radius is taken to be 
equal to one half of the half-mass radius of the system. The symmetrization is performed by 
accepting the radius and the magnitude of the velocity of each simulation particle, following 
the procedure to generate clumpy conditions discussed in TBvA05, and by redistributing 
uniformly the angular variables in both position and velocity space. This initialization pro- 
cedure leads to a smooth initial density profile, decreasing approximately linearly in radius 
(p(0) ~ 2p{rM))j which creates a potential well that is deeper than the one of a uniform 
sphere with the same cutoff radius. The initial system is non-rotating and isotropic. Mild 
correlations in the magnitude of the velocities are present, since there is a residual memory 
of the cold clumpy state. The precise form of the density profile as well as the strength of 
the velocity correlations depend on the details of the initial random positions of the clump 
centers, i.e. on the seed of the random number generator. 

A snapshot of typical initial conditions from which our set of simulations of coUisionless 
collapse starts is reported in Fig. 1. The correlations in the magnitude of the velocities, char- 
acteristic of the original clumps from which the present initial conditions are extracted, can 
be seen as ring patterns created by the particles in the {vx, Vy) space. We should emphasize 
that these initial conditions are rather artificial and have been employed with the specific 
purpose of studying the effects of the suppression of mixing in phase space. It appears un- 
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likely that in a realistic scenario of galaxy formation the initial state is characterized by such 
high degree of symmetry. Nevertheless, the properties of the end-products of these collapse 
simulations are not too unrealistic, at least in their outer parts (see Subsection 2.1). 

In Table 1 we summarize the properties of this set of simulations. The initial conditions 
of the SI series are generated using the same positions for the clump centers but a different 
total number of particles (to ensure that the properties of the end products do not depend 
on A^). Run 5*1^ has 8x 10^ simulation particles; 5*1* and Sl^ have 10^ particles, but the first 
simulation is run with our particle-mesh code, while the second with Dehnen's code. Runs 
Sa and Sb arc colder versions of 5*1* (obtained by means of a global rescaling of the initial 
velocities by a constant factor). Runs S2 and S3 are generated using a different seed for the 
random numbers. In order to characterize the level of anisotropy achieved in the end-states 
obtained from the simulations, we refer to the global anisotropy parameter 2Kr/ Kt and to 
the anisotropy profile a{r), defined as a{r) — 2 — {{vq) + {v'^))/{v^). 

During collapse, the spherical symmetry is well preserved, as can be seen for SI not only 
from the evolution of the eigenvalues of the inertia tensor of the system (Fig. 2), but also 
from the conservation of the single particle angular momenta (Fig. 3). Mass loss (i.e. the 
number of particles that acquire a positive energy during the collapse) is limited, well below 
the loss recorded for homogeneous uniform spheres with similar initial virial ratio u, where 
the system can lose up to one third of its total mass. The combination of spherical symmetry 
and limited mass loss leads to high final density concentrations, with p(0)/p(rM) ^ 1500 in 
run SI. As shown in Fig. 4, the density profile is reasonably well represented by a rather 
concentrated f^'^^ model (see TB05) or by a Jaffe density profile (JafFe 1983). 

The global amount of pressure anisotropy (see Fig. 2) evolves rapidly in the first few 
dynamical times and then reaches its quasi-equilibrium value at t ^ St^^. When 2Kj./Kt 
reaches its peak value, at the time of maximum contraction of the system, the anisotropy 
radius (defined implicitly by a{r-a) = 1) is located well inside: at that time, for run SI 
the mass within Tq is only 1% of the total mass, while later the sphere associated with the 
anisotropy radius contains approximately 20% of the total mass. The final global content 
of pressure anisotropy is high also for run 5*3, which starts from moderately warm initial 
conditions {u — 0.25; for comparison, see the results of non-symmetrized runs with phase- 
space mixing presented in TBvAOS). 

Interestingly, the central regions have a final pressure anisotropy profile slightly biased 
toward tangential orbits (see Fig. 4). This effect appears in the high resolution <S'l''" simu- 
lation, with 8 X 10^ particles; the realization with 10^ particles (5*1*) does not exhibit this 
feature and indeed is characterized by a slightly higher value of 2Kr/KT- In any case the 
transition from isotropic to radial pressure is very sharp. The shape of the anisotropy pro- 
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file cannot be represented either by a profile similar to those of the /^^^^ models or by an 
Osikpov-Merritt profile (Osipkov 1979; Merritt 1985). 

To check the robustness of the numerical results, with Dehnen's tree code we have let the 
final configuration reached in 5*1* evolve for 30 additional dynamical times, without noticing 
any sign of significant changes. In addition, we have run a collapse simulation starting with 
the same initial conditions as 5*1* using Dehnen's tree code during the entire simulation 
(nm 5*1^); this simulation shows no macroscopic differences with respect to 5*1* (these are 
summarized in Table 1 and are at the level of 1% ). 

If we consider colder and colder initial conditions within the framework of simulations 
considered in this paper, the radial orbit instability eventually sets in. A reduction of the 
initial virial ratio u below 0.05 for Sl-like initial conditions (by rescaling the velocities 
by a constant factor) leads to runs that show evidence for the radial orbit instability: a 
simulation with u = 0.03 {Sb) leads to a strongly flattened system, with a final amount of 
global anisotropy close to 2.5. 

A simulation starting from u = 0.05 shows an interesting behavior which we may identify 
as that of marginal stability, since the final state is characterized by an aspect ratio (see 
Fig. 5) 7] that oscillates between 0.99 and 0.93, on a time scale longer than the dynamical 
time; the related anisotropy content is high {2Kr/KT pa 2.75). 

2.1. Elliptical galsLxies and the properties of the end-products of some 
symmetric colUsionless collapse simulations 

As we have seen above, the end-products of this set of simulations, starting from highly 
symmetric initial conditions, are characterized by a density profile which is not too far 
from that of models, such as the Jaffe (1983) profile, that are known to be associated with 
realistic surface brightness profiles for bright elliptical galaxies (under the assumption of 
constant mass-to-hght ratio). To quantify the effects of the systematic differences between 
these models and the products of the simulations, in Fig. 6 we consider the projected density 
profile measured at the end of simulation Sl~^ and we compare it to the i?^/^ (de Vaucouleurs 
1948) and to the i?-*^/" (Sersic 1968) laws. The residuals from the best fit, at fixed effective 
radius Rg, are not too large, especially for the Sersic law with n ^ 8.3. In fact, the deviations 
are within 0.2 magnitudes over a wide interval in radius, from 0.2 to 10 Rg, however, in the 
central regions the difference can be as high as 0.5 magnitudes, probably because of the 
presence of the "bump" in the density profile around 0.1 tm (see Fig. 4). Residuals from the 
i?^/^ are larger. 
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For a direct comparison of the models considered in this paper and the observations, one 
might think of making a kinematical test. Unfortunately, as to the anisotropy of the velocity 
dispersion tensor, at present a comparison with the observations would be extremely difficult 
to carry out. The best empirical evidence for radial anisotropy in elliptical galaxies comes 
from the observed flattening (and triaxiality) in objects that are not rotationally supported. 
In principle, pressure anisotropy affects the shape of the velocity distribution integrated 
along the line of sight, so that its presence could be detected by observing the shapes of 
the line profiles (and in particular the deviations from Gaussian profiles). In practice, the 
effects on line profiles are rather modest even for significant amounts of pressure anisotropy, 
so that, within the uncertainty limits of the measurement, the pressure anisotropy profile of 
elliptical galaxies remains basically undetermined by current observations (e.g., see Gerhard 
1993; Gerhard et al. 2001; dc Zccuw et al. 2002). In other words, it would be very hard to 
ascertain whether some ellipticals are indeed as isotropic in their inner regions as implied by 
some of the models that wc have investigated. 

As briefly mentioned in the Introduction, numerical simulations of galaxy formation 
suggest that elliptical galaxies may be close to marginal stability with respect to the radial 
orbit instability, close to the threshold stated by Polyachenko & Shukhman (1981). In this 
respect, studying the effects on the radial orbit stability of a highly isotropic core, as we do 
in this paper, is of interest, because there are astrophysical processes that are often ignored 
in the simulations of galaxy formation and should operate in real elliptical galaxies, which 
may favor a more isotropic core or even an anisotropic core with a slight bias of the pressure 
tensor in the tangential directions. Here we have in mind processes such as those due to 
the presence of a super-massive black hole (for the effects of a black hole on the pressure 
anisotropy within the sphere of influence see e.g., CipoUina & Bertin 1994, Baumgardt et al. 
2004). To be sure, for a physically signiflcant application of these ideas, one should check 
the time-scales of the competing processes that are invoked, before drawing the relevant 
conclusions. 



3. Exceptionally stable equilibria constructed from the Jeans equations 

We now address the issue of whether a given initial conflguration, characterized by 
assigned density p(r) and pressure anisotropy a{r) proflles, is stable with respect to the 
radial orbit instability. Wc arc not aware of models with analytical distribution function able 
to incorporate the sharp feature in the anisotropy profile, of the kind observed at the end 
of the simulations described in the previous Section. Therefore, we decided to initialize the 
simulations by means of candidate equilibrium solutions obtained from the Jeans equations. 
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as outlined in the Appendix. We should emphasize that the simulations that we describe 
below in this Section are simulations of candidate equilibria. In fact, since we never reach 
the point of actually reconstructing an underlying equilibrium distribution function, if we 
happened to find significant evolution we could be either in a situation of genuine instability, 
or, more simply, in a situation of non-equilibrium. In turn, since we will show cases where 
wc do not find such evolution, we may claim that indeed we have found not only a genuine 
quasi- equilibrium state but also proved that it is approximately stable. 

To model a density profile of the kind found at the end of the 5*1 simulations, when no 
evolution is observed over several dynamical times, we use a superposition of a regularized 
Jaffe (1983) profile: 

with A, e, and b free scales, and a central core of the form: 



Pci^) = (^2 + ^^2)20 ' (2) 

where A' is a free scale, while ^ is a dimensionless parameter of order 1 ~ 0.6). The form 
for the density Pmod = Pj + Pc is taken for convenience, so as to reproduce not only the 
large-scale structure of the density distribution realized in the SI simulations, but also the 
bump in the density profile around tm/S (see Fig. 4). 

As noted earlier, the SI anisotropy profile is fiatter than the corresponding Osipkov- 
Merritt profile with same at small radii, while it is less steep at large radii. Thus we chose 
to represent the profile with: 



a{r) = 2 ^ (3) 

with Ta and 7 being free parameters. A single choice of the (r^; 7) values is unable to correctly 
reproduce the anisotropy profile measured from the simulation over the entire radial range. 

Thus wc fit separately the anisotropy in the "core", up to a radius Vch ^ ^"a, where 7 = 4, 
and in the "halo", i.e. for r > Vch, where 7 = 6/5. In a neighborhood around Vch the two 
profiles are matched so that the final a{r) and its derivative are continuous functions. 

With a suitable choice for the various parameters that define the above functions, the 
profiles obtained from the simulation can be fitted with an accuracy of better than 10 percent. 
We take this as a good starting point to investigate the stability of equilibrium configurations 
similar to those produced in the simulations of coUisionless collapse. 
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We have first studied the evolution of a model initialized with density and anisotropy 
profiles similar to those of runs 5*1 (see entry Jl in Table 2 and Fig. 7). We have then 
proceeded to study the evolution of neighboring configurations by slightly modifying the 
density and/or the anisotropy profiles. In particular, we have considered models with a den- 
sity profile without the inner "bump" (i.e., without the pc contribution), and different forms 
of anisotropy profile, ranging from steep profiles over the entire radial range, to Osipkov- 
Merritt and /(''^-hke profiles. Interestingly we have found that although none of the various 
combinations turns out to be violently unstable, nevertheless, the only simulation where 
practically no sign of evolution occurs is Jl, the one associated with the profiles that best 
fit those of 5*1 . 

For the systems that show definite signs of evolution, the evolution appears to be very 
slow. For example, in the J2 simulation the initial configuration lasts basically unchanged 
for more than 10 dynamical times. In this case, illustrated in Fig. 8, the system remains 
close to spherical symmetry with 2Kr/KT ~ 2.9, and ends up only much later as a prolate 
system. 

4. Conclusions 

In this paper we have studied and clarified the dependence of the radial orbit instability 
on the shape of the anisotropy profile in the central region of a stellar system. By means 
of a scries of numerical simulations we have shown that stable, centrally isotropic equilibria 
with a significant global amount of anisotropy can be reached during highly symmetric cold 
collapse events or initialized by solving the Jeans equations. In particular we have found a 
metastable state with 2Kr/ Kt ~ 2.9. 

The experiments that wc have performed suggest that the presence of an isotropic 
core may act as an important stabilizing factor for the radial orbit instability. The two- 
component (core-halo) structure for run Sl^ is indeed evident not only in the density and 
pressure anisotropy profiles (see the structure out to r ^ in Fig. 4), but also in the 

phase space distribution N{E), as indicated by the peak located at high values of the binding 
energy {E fa —14, in code units; see left panel of Fig. 9). 

One indication for the stabilizing effect of a highly isotropic core may be found by 
recalling that the linear modal stability analyses have shown that the density profiles of the 
fastest growing modes are peaked in the innermost region of the system (see e.g., Bertin 
et al. 1994, Fig 8). Thus we may argue that, if the core were almost perfectly isotropic, 
there would be no room for the modes to be excited. 
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These results are qualitatively reminiscent (although the physical mechanisms that are 
involved are completely different), of the stabilizing role that a hot bulge can provide in 
relation to the stability of sclf-gravitating disks, as noted in a number of papers starting 
with Berman & Mark (1977) and Scllwood (1981). For disks, the detailed mechanisms 
underlying the origin or the suppression of global bar and spiral instabilities are reasonably 
well understood and known to depend on the key structural properties of the basic state 
(i.e., the disk density, effective velocity dispersion, and differential rotation profiles). For 
anisotropic spherical stellar systems we still lack a clear picture of the relevant underlying 
mechanisms. In this respect, this paper offers an interesting clue to more systematic studies 
that should be devoted to investigating the nature of the radial orbit instability in two- 
component systems, also in view of the central properties of the dark halo. 

We would like to thank T.S. van Albada for a number of conversations and useful 
suggestions. This work has been partially supported by the Italian Ministry MIUR. 

A. Candidate collisionless equilibria generated by the Jeans equations 

The goal is to construct stellar dynamical configurations, as initial conditions for N-body 
simulations, corresponding to systems with given density and pressure anisotropy profiles. 
As is well known, in general this problem is not well-posed, because a given pair (p(r), a{r)) 
may even lack a physically acceptable underlying equilibrium distribution function. We 
proceed by applying the necessary constraints of the Jeans equations and then initialize our 
simulations with a procedure described below, well aware that we may actually miss the 
desired equilibrium conditions (see comment at the beginning of Sect. 3). 

The knowledge of the density profile immediately identifies the self-consistent potential 
$(r). In addition, it defines the integrated mass profile, from which the positions of the 
particles can be correctly initialized. In the absence of an analytical distribution function, to 
complete the initialization of the N-body collisionless candidate equilibrium configuration we 
then resort to the Jeans equations to extract the information about the velocity dispersion 
profile that would be required by the conditions of equilibrium. In practice, for a non-rotating 
spherically symmetric system the relevant hydrostatic equilibrium equation can be written 
as: 



where $ represents the mean field gravitational potential (which can be taken to be known 




(Al) 
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for known p(r)) and cr^ = {v^) the velocity dispersion in the radial direction. For assigned 
profiles (p(r), Q;(r)), under the natural boundary condition given by 

p(oo)(7,2(oo) = 0, (A2) 

the equation can be solved for cr^: 

Thus we have obtained the velocity dispersion profiles consistent with the given density and 
anisotropy profiles. 

At this point, we initialize our N-body code by assuming that the velocity dispersion 
profiles obtained above correspond, at least approximately, to a truncated Gaussian distri- 
bution in velocity space (the truncation is imposed by the requirement that only bound 
particles, i.e. particles with negative energy, belong to the system). This assumption gives a 
reasonable description of systems for which the velocity distribution is determined by physi- 
cally motivated distribution functions (e.g., for the stable f^^"* models); but, in general, there 
is no guarantee that N-body systems initialized by this method be in approximate dynamical 
equilibrium. 

In our case, the limitations of this approach are not severe, because we are interested in a 
stability study. If we happen to obtain a configuration with no significant signs of evolution, 
we are certain a posteriori that a quasi-equilibrium state has indeed been produced. On the 
other hand, even if the initialization gives a configuration not in dynamical equilibrium, but 
a nearby equilibrium configuration is approached on a very short time scale (typically on 
the order of one dynamical time), then we may argue that we have also obtained our goal 
of producing a quasi-equilibrium state with characteristics close to those desired. 

A.l. Tests on known distribution functions 

Wc have tested our method on isotropic systems by initializing Plummer and polytropic 
models and then on some anisotropic configurations associated with the /'-'^^ family (see 
TB05), under conditions of stabihty. By construction, we have a good match at the level of 
density, velocity dispersion, and anisotropy profiles. Differences in the fourth order moment 
of the velocity distribution are typically at the level of a few percent. The models initialized 
with the method outlined above appear to be stable, except for some modest re- arrangements 
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of the initial anisotropy profile in the case of the /^'^^ models. In the latter case, there is a 
tendency for the anisotropy content in the central region to decrease so that a flatter profile 
in a{r) is eventually generated. 
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Fig. 1. — Snapshot of the initial conditions in the planes {x, y) (left panel) and {vxi Vy) (right 
panel) for simulation SI. The plots are given in code units. The effects of the symmetrization 
over the initial conditions of the kind considered in TBvAOS are easily appreciated in velocity 
space, where dumpiness before symmetrization was high (see Fig. 3 in TBvAOS); here 
particles tend to form ring structures. The adopted symmetrization of initial conditions 
suppresses mixing in phase space during coUisionless collapse. 
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Fig. 2. — Left panel: evolution of the ellipticities e — b/a and r] = c/a, where a>b>c are 
the values of the axes of the system evaluated from the inertia tensor, for the Sl'^ simulation. 
Right panel: evolution of the ratio 2Kr/ Kt for the same simulation. The dynamical time is 
defined in terms of the total mass and of the total energy ostd — GM^/"^ / {—2Etotf^'^ ■ 
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Fig. 3. — Left panel: variation of one component of the specific angular momentum of the 
simulation particles for run Sl~^. Since during evolution the system remains spherically 
symmetric, the single particle angular momenta are approximately conserved. Right panel: 
Distribution of initial angular momenta (solid line) for the same simulation: the distribution 
basically coincides with that for the final angular momenta (dotted line). 
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Fig. 4. — Density (in code units; left) and anisotropy profile (right) for the final state reached 
in the Sl~^ simulation. The solid line in the left panel is the density profile associated with 
the (1/2; 6.4) f^'^^ model (see TB05 for definitions and notation). The anisotropy profile 
is compared to the Osipkov-Merritt anisotropy profile (solid line) and to the profile for an 
model characterized by a similar amount of global anisotropy (dashed line; the (1/2; 1) 
model, which has the same as the end-state for ^l"*"). Clearly the theoretical models are 
unable to capture the rapid increase in a at r 
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Fig. 5. — Ellipticities (left) and anisotropy content (right) for the Sa simulation, as plotted 
in Fig. 2. Note the extremely high value of anisotropy content {2Kr/KT > 20) reached 
during the collapse. 



- 21 - 



0.5 



0.5 




0.5 



1 



1.5 



R/R 



Fig. 6. — Residuals from the fit with the R^^^ law (sohd line) and with the i?^/" law (dashed 
line), in magnitudes, of the projected density profile for the final configuration reached in 
simulation Sl'^ . The fit has been obtained by considering a constant mass-to-light ratio, at 
fixed effective radius Rf.- The best fitting index n \s n ^ 8.3. Note the sizable residuals in 
the central regions. 
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Fig. 7. — EUipticities (left) and anisotropy content (right) for the Jl simulation, as plotted 
in Fig. 2. 
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Fig. 8. — EUipticities (left) and anisotropy content (right) for the J2 simulation (2X^1 Kt — 
2.92 at i = 0), as plotted in Fig. 2. 
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Fig. 9. — Phase space distributions N{E) (left) and N{E, .P) (right) for the final state 
reached in the 5*1+ simulation. The peak in the bottom left panel around E —14 (in code 
units) is associated with the isotropic core. 
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Table 1: Collapse simulations. The simulations have been run with N — 10^ particles, except 
5"!+ for which N ^ 8 x 10^. The columns list the initial virial ratio u = {2K/\W\)t=Q, the 
fractional mass loss AM = [M{t = 0) — M{tend)]/M{t — 0), the final global anisotropy 
K. = 2Kr/ Kt of the system of bound particles, its anisotropy radius (i.e. the radius where 
the local anisotropy parameter is a = 1) in units of the half-mass radius, and the cllipticities 
t = h/a and 7] = c/a, where a ^ b ^ c are the axes of the final quasi-equilibrium state 
computed from the inertia tensor. 
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Table 2: Simulation runs initialized with models constructed from the Jeans equations. The 
simulations have 2 x 10^ particles. Here we list the initial and final global anisotropy k = 
2Kr/KT of the system, its final anisotropy radius in terms of the half-mass radius, and the 
final ellipticities e — h/a and r) = c/a, where c are the axes evaluated from the inertia 

tensor. The initial conditions are summarized in the last column and range from the best 
fit of 5"!, simulation Jl, to regularized Jaffe profiles pj (see Eq. 1) with Osipkov-Merritt or 
f^^^ like anisotropy. Note that simulation J4, initialized with a regularized Jaffe density plus 
an unstable f^^"* anisotropy profile {k ~ 2.3) evolves rapidly within the first dynamical time 
and the anisotropy is quickly reduced below k — 2 while preserving the spherical symmetry. 
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